International Journal of Modern Physics C, 
© World Scientific Publishing Company 



IMPROVED SPIN DYNAMICS SIMULATIONS 
OF MAGNETIC EXCITATIONS 



D. P. LANDAU, SHAN-HO TSAI, M. KRECH, and ALEX BUNKER* 
Center for Simulational Physics, The University of Georgia, 
Athens, GA 30602, U.S.A. 



Received (received date) 
Revised 

Using Suzuki- Trotter decompositions of exponential operators we describe new algo- 
rithms for the numerical integration of the equations of motion for classical spin sys- 
tems. These techniques conserve spin length exactly and, in special cases, also conserve 
the energy and maintain time reversibility. We investigate integration schemes of up to 
eighth order and show that these new algorithms can be used with much larger time 
steps than a well established predictor-corrector method. These methods may lead to 
a substantial speedup of spin dynamics simulations, however, the choice of which order 
method to use is not always straightforward. 

1. Introduction 

Our understanding of static behavior near phase transitions is now mature and has 
resulted largely from the investigation of simple model spin systems such as the 
Ising, the XY, and the Heisenberg model. These models are equally valuable for 
the investigation of dynamic critical behavior and dynamic scaling. Realistic models 
of magnetic materials can be constructed from these simple spin models, however, 
the theoretical analysis of experimentally accessible quantities, such as the dynamic 
structure factor, is usually too demanding for analytical methods. Computer sim- 
ulations are beginning to provide important information about dynamic critical 
behavior and material properties of model magnetic systems Ercrtil. These simula- 
tions use model Hamiltonians with continuous degrees of freedom represented by a 
three-component spin S& with fixed length |S&| = 1 for each lattice site k. A typical 
model Hamiltonian is then given by 

H = -J £ (S%Sf + S%S? + ~Dj2(S z k f , (1) 

<*:,(> k 

where J is the exchange integral, < k,l > denotes a nearest-neighbor pair of spins 
Sfc, A is an anisotropy parameter, and D determines the strength of a single-site or 
crystal field anisotropy. 

Modelling specific magnetic materials may require additional interactions in the 
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Hamjltonian, such as two-spin exchange interactions between more distant neigh- 
bors u, three spin exchange interactions, or even biquadratic coupling 

The thermodynamic properties can be obtained from a Monte-Carlo simulation 



and the dynamical properties, 
equations of motion given by 



he spin system are provided by solutions to the 



These equations must be integrated numerically, where a Monte-Carlo simulation of 
the model provides equilibrium configurations as initial conditions for Eq.(^). The 
most important quantity to be extracted from the numerical results is the dynamic 
structure factor S(q, oS), which is given by the space-time Fourier transform of the 
spin-spin correlation function 

g^{v k -T h t-t') = {s^t)s?{t% (3) 

where a, (3 — x,y, z, rj, and r/ are lattice vectors, and the average (. . .) must be 
taken over a large number of independent initial equilibrium configurations. (This 
procedure is appropriate since the typical time over which Eq. (|J) can be integrated 
is much shorter than typical timescales set by other excitations.) 

To speed up the numerical integration of Eq. (||) it is desirable to use the largest 
possible time step; however, with standard methods the size of the time step is 
severely limited by the accuracy within which the conservation laws of the dynamics 
are obeyed. It is evident from Eq.(|^) that the total energy is conserved, and if, 
for example, D = and A = 1 (isotropic Heisenberg model) the magnetization 
M = J^k Sfc is a l so conserved. For the anisotropic Heisenberg model, i.e., A ^ 1 or 
D =/= only M z is conserved. Conservation of spin length and energy is particularly 
crucial, because the condition |S^ | = 1 is a major part of the definition of the model 
and the energy of a configuration determines its statistical weight. It would therefore 
also be desirable to devise an algorithm which conserves these two quantities exactly. 

In the remaining sections we describe a 4th-order predictor-corrector method 
and a new integration procedure, which is based on Trotter-Suzuki decompositions 
of exponential operators Bfl. We compare both schemes with special regard to speed 
and the accuracy within which the conservation laws hold. 

2. Integration Methods 



2.1. Predictor- corrector methods 

Predictor-corrector methods have been quite effective for the numerical integra- 
tion of spin equations of motion; however, in order to limit truncation errors small 
time steps St must be used with at least a fourth-order scheme. In a more symbolic 
form Eq.(|^) can be written as y = f(y) with the initial condition y(0) — yo, where 
y is a short-hand notation of a complete spin configuration. The initial equilibrium 
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configuration is denoted by yo- The predictor step of the scheme used here, the 
explicit Adams-Bashforth four-step method H, is 

y(t + St) = y(t) + ^ [55/(y(f)) - 59/(y(i - St)) + 37f(y(t - 2St)) - 9f(y(t - 3St))} 

(4) 

which has a local truncation error of the order (St) 5 . The corrector step consists of 
typically one iteration of the implicit Adams-Moulton three-step method El 

y(t + 5t) = V(t) + Y A [9/(y(i + Si)) + W(y(t)) ~ 5f(y(t - St)) + f(y(t - 2St))} (5) 

which also has a local truncation error of the order (St) 5 . (Values for y(5t), y(2St), 
and y(35t) in addition to y(0) = yo can be provided by three successive integrations 
of y = f(y) (see Eq.(|^)) by the 4th-order Runge-Kutta method □ .) This method 
requires that spin configurations at the last four time steps must be kept in memory. 

This predictor-corrector method is very general and is independent of the spe- 
cial structure of the right-hand side of the equations of motion (see Eq.(|^)). The 
conservation laws discussed earlier will only be observed within the accuracy set by 
the truncation error of the method. In practice, this limits the time step to typically 
St = 0.01/ J in d = 3 for the isotropic model (D — 0) 13, where the total integration 
time is typically 600/ J or less. 

2.2. Suzuki- Trotter decomposition methods 

The motion of a spin may be viewed as a precession of the spin S around an 
effective axis il which is itself time dependent. The lattice can be decomposed into 
two sublattices such that a spin on one sublattice precesses in a local field £1 of 
neighbor spins which are all located on the other sublattice. For the Hamiltonian in 
Eq.(Q) there are only two such sublattices if the underlying lattice is simple cubic. 

To illustrate the method, we consider first the D = case. The basic idea of the 
algorithm is to rotate a spin about its local field f2 by an angle a = \Q\5t, rather than 
directly integrate Eq.(||). This procedure guarantees the conservation of the spin 
length |S| and energy to within machine accuracy. Denoting the two sublattices by 
A and B, respectively, we can express the local fields acting on the spins on sublattice 
A and B as f2^i[{S}] and £!g[{S}], respectively. The set of equations of motion for 
spins on one sublattice reduces to a linear system of differential equations if the 
spins on the other sublattice are kept fixed. Thus an alternating update scheme 
may be used, i.e., the spins Sk^A are rotated for the given values of S^gg and vice 
versa. (The scalar products S/cg^ ■ f2_4[{S}] remain constant during the update of 
Sfc e .4 and the scalar products SkeB ' remain constant during the update of 

Sfcgg). Note, that each sublattice rotation is performed with the actual values of 
the spins on the other sublattice, so that only a single copy of the spin configuration 
is kept in memory at any time. However, the magnetization will not be conserved 
during the above rotation operations. Since the two alternating rotation operations 
do not commute, a closer examination of the sublattice decomposition of the spin 
rotation is required. 
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We now decompose a full configuration into two sublattice components y_\ and 
yg, i.e. y — (y_4, yg), and denote by matrices A and -B the generators of the rotation 
of the spin configuration y.4 on sublattice A at fixed yg and of the spin configuration 
yg on sublattice B at fixed y_4, respectively. The update of the configuration y from 
time t to £ + St is then given by an exponential (matrix) operator 

y(t + &) = e {A+B ^ t y{t). (6) 

Although the exponential operator in Eq.(^) rotates each spin of the configuration 
it has no simple explicit form, because the rotation axis for each spin depends on the 
configuration itself; however, the operators e A5t and e BSt which rotate y.4 at fixed 
yg and yg at fixed y.4, respectively, do have a simple explicit form. We demonstrate 
this for the case A = 1 and D = in Eq.(|l|). For each keAwe find 

^[{s}] = -j s < = ^> ( 7 ) 

where NN(k) denotes the nearest neighbors of k (which belong to yg). Eq.(^) can 
be readily generalized for A ^ 1, the case D ^ will be discussed below. The 
explicit rotation of spins on each sublattice reads (see also Watson et al a) 



n k (n k -s k (t)) 

S k (t+6t) = -2 



^k(t) ^2 



cos(|S2fc|d£)H tt—. sin(|S2 fe 



(8) 

Note that according to Eq.(||) f2fc • Sfc(f + 5i) = fl k ■ S k (t) and energy is thus 
conserved. The alternating update scheme for the integration of the equations of 
motion amounts to the replacement e ( A + B ) St — > e -4<5t e _Bi5t j n Eq.(§), which is only 
correct El up to order (St) 2 . The magnetization will therefore only be conserved up 
to terms of the order St (global truncation error), but one can employ higher order 
Suzuki- Trotter decompositions of the exponential operator in Eq.(^) to decrease 
the local truncation error of the algorithm and thus improve the conservation. The 
simplest possible improvement is given by the 2nd-order decompositior£2l 

e (A+B)8t = e A8t/2 e B8t e ASt/2 + ^3) (g) 

which is equivalent to the midnoint integration method applied to the equations of 
motion (see also Watson et al u) . We can also use the mth-order decomposition E3 

n 

e (A+B)St = -Q ePl A8t/2 ePl B8t ePl A8t/2 + Q^m+l) (1Q) 
i=l 

where n = 5 for 4th-order and n — 15 for 8th-order and the parameters p, are given 
by Suzuki and Umeno 113. 

The additional computational effort needed to evaluate higher order expressions 
can be compensated to some extent by using larger time steps. The evaluation of the 
trigonometric functions in Eq. (B) can also be avoided since the above decompositions 
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are only correct to within a certain order in St and it is therefore sufficient to replace 
sin a; and cosx by appropriate Taylor polynomials; alternatively a Cayley transform 
could be used (up to 2nd and 4th-order, it corresponds to sin a; = x p(x)/[p 2 (x) + 
(x/2) 2 ], with p(x) — 1 and p(x) = 1 — a; 2 /12, respectively; determining cosx from 
sin 2 x + cos 2 x = 1 ensures spin- length conservation) . For a 4th-order method the 
Cayley transform was 10 — 20% faster, depending upon the machine. Note that the 
decompositions maintain the time inversion property of e ( A + B ) St . The inclusion of 
next-nearest neighbor bilinear interactions on a simple cubic lattice can be treated 
within the above framework if the lattice is decomposed into four sublattices. 

This approach can also be extended to the case D ^ 0, but in contrast to the 
isotropic case, the equation of motion for each individual spin on each sublattice is 
nonlinear. In practice, the best form of solution is via iterative numerical methods. 
For the sublattice decomposition of the spin rotation the requirement for energy 
conservation in the presence of a single site anisotropy is 

n k ■ S k (t + St)-D [S z k (t + St)} 2 = fi fc • S fc (t) - D [S z k (t)f (11) 

for k € A and k 6 B, where fl k is given by Eq.(Q). In order to perform a rotation 
operation in analogy to Eq.(^) we have to identify an effective rotation axis. This 
can be achieved by rewriting Eq. ([il]) in the form fl k ■ (S k {t + St) — Sfc(t)) = 0, where 

n k = n k -D (o, o, si® + s z k (t + st)) . (12) 

Since the rotation requires knowledge of S k at the future time t+St, this problem can 
be solved iteratively starting from the initial value SI (t+8t) = S k (t)+(U k xSk(t)) z St 
in Eq.(p"2) and performing several updates according to the decompositions given 
by Eqs^(| ) or (|To|), respectively, in order to improve energy conservation according 
to Eq.rtTj). Both the degree of conservation and the execution time depend to some 
extent on the number of iterations used. The initial value for SI (t + St) used here 
yields a better energy conservation (at almost no extra CPU time) than using the 
initial value Sl(t + St) = S k (t), with the same number of iterations. Biquadratic 
interactions can be treated by the same iterative scheme, but inclusion of three-spin 
interaction would require reconsideration of the sublattice decomposition. 

3. Results and Comparisons 

For a quantitative analysis of the integration methods outlined above we restrict 
ourselves to the Hamiltonian given by Eq.(|l|) for A = 1 in d = 3. The underlying 
lattice is simple cubic with L — 10 lattice sites in each direction and periodic 
boundary conditions in all cases discussed below. 

In order to compare the different integration methods we first investigate the 
accuracy within which the conservation laws are fulfilled. The initial configuration 
is a well equilibrated one from a Monte-Carlo simulation for A = 1 at a temperature 
T = 0.8T C for D = and D = J, where T c refers to the critical temperature 
of the isotropic model (D = 0). The magnetization of such a configuration is 
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non-zero and provides an indicator for the numerical quality of the magnetization 
conservation. We integrate the equations of motion to t — 800/ J and monitor the 
energy e(t) = E(t)/ (JL 3 ) of the configuration per spin and the modulus m(t) = 
\M.(t)\/L i of the magnetization per spin for the isotropic case D = and its z- 
component m z (t) = M z (t) / L 3 for the strongly anisotropic case D = J as functions 
of time. Note that for these tests both integration methods are started from identical 
initial configurations. 

For D = the implementation of Eqs.([}j) and (|Io| ) using Eq.(^) is straightfor- 
ward. The Taylor polynomial for sin a; is chosen as sin a; = x — x 3 /6 for Eq.(^) 
and up to (and including) the terms x 5 /120 and x 9 /9! for the 4th- and 8th-order 
decompositions in Eq.(jl0|), respectively, in order to reduce the magnetization fluc- 
tuations. Fig.[j] shows that e(t) for the predictor-corrector method increases lin- 
early with time whereas the decomposition methods both yield e(t) — const. Thus, 
St = 0.01/ J is about the largest value that can be used without introducing substan- 
tial non-conservation of the energy. Fig.|] displays the magnetization conservation 
for the 2nd-, 4th- and 8th-order decomposition methods, all with the same time 
step St = 0.1/ J. The predictor-corrector method conserves m(t) exactly, whereas 
the decomposition methods cause fluctuations of m(t) on all time scales. It is also 
clear that the second-order method is unstable with such a large time step. The 
temporal structure of m(t) for the eighth-order decomposition methods with differ- 
ent time steps is displayed in Fig.||. Even with a time step as large as St = 0.25/ J, 
the magnetization is rather well conserved out to a time of t max — 800J -1 . The 
three different decomposition methods can be made to produce fluctuations of the 
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Fig. 1. Energy e(t) = E(t)/{JL 3 ) per spin for the predictor-corrector method (see Eqs.(^) and 
(^)) for several time steps and D = 0. 

same magnitude, e.g. of 2 x 10 -5 , by adjusting the value of St to be 0.007J -1 for 
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Fig. 2. Magnetization m(f) = |M(t)|/L 3 per spin for different order decomposition schemes for 
D = and time step St = 0.1/ J: (dotted line) 2nd-order scheme; (solid line) 4th-order scheme; 
(dashed line) 8th-order method. 
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Fig- JL, Magnetization m(t) per spin for D = using the eighth-order decomposition method (see 
Eq.(hoj)) with different time steps. 
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the 2nd-order method, 0.1J -1 for 4th-order and 0.25J -1 for 8th-order. A single 
integration of the equations of motion using Eq.(^|) (2nd-order decomposition) is 
about twice as fast as the predictor-corrector method. The 4th- and 8th-order de- 
compositions (see Eq.(|lO|)), however, are respectively about 2.5 and 9 times slower 
than the predictor-corrector method. Taking the change in time step by factors 
of 0.7, 10, and 25, respectively, into account, we find that the 2nd-, 4th- and 8th- 
order decomposition methods yield a speedup of the integration of the equation of 
motion by factors of approximately 1.5, 4, and 2.5, respectively. Using time steps 
St = 0.04/J and 0.2/ J for the 2nd- and 4th-order methods, respectively, yields in 
both cases an eightfold speedup as compared to the predictor-corrector method with 
St = 0.01/ J. If the overall quality of the magnetization conservation is also taken 
into account, there is a clear advantage for the 4th-order decomposition according 
to Eq.(|ic|) for the isotropic case D = 0. Somewhat surprisingly, with the current 
implementation, the 8th-order method is not competitive. For strong anisotropy, 
D = J, the predictor-corrector method can be applied as before, but the decompo- 
sition scheme must be modified because the spin rotation axis depends on the spin 
value S| at the future time t+5t (see Eq.([l2"|)). As described in Sec. 2.2 this self con- 
sistency problem is solved iteratively, where the quality of the energy conservation 
depends on the number of iterations performed. For the 2nd-order decomposition 
with St = 0.04/J two iterations are sufficient to obtain a better energy conserva- 
tion than the predictor-corrector method. Thus a time integration step using the 
2nd-order decomposition takes twice as long as for D = 0, so that its advantage in 
speed only comes from the increase in the time step, which is still a factor of four. 
For the 4th-order decomposition with St = 0.2/ J six iterations are needed to obtain 
energy conservation to within six significant digits, so one integration becomes 15 
times slower than one with the predictor-corrector method. From the increase of 
St by a factor 20 only a 30% gain in speed is obtained. The number of iterations 
needed decreases with St, but this decrease does not compensate the loss in speed 
due to the smaller time step. However, one still obtains a greatly improved energy 
conservation. All methods behave similarly, the change in energy is basically linear 
with time (see Fig.||) . The reason for this is the iterative nature of all four methods 
in the case D ^ 0. The overall accuracy of the magnetization conservation appears 
to be independent of D for all decomposition methods. Considering both overall 
energy conservation and speed, the 2nd-order decomposition has some advantages 
over the predictor-corrector method. 

Lastly, in Fig.|5| we show spin dynamics data for the dynamic structure factor in 
the (100) direction for the isotropic Heisenberg ferromagnet (D = 0) on a simple 
cubic lattice (L = 10) at T = 0.8T C obtained with the second-order decomposition 
method (see Eq.([|), St — 0.04/J). The equations of motion have been integrated 
to 800/ J and averages have been taken over 1000 initial configurations, where the 
time displaced correlation functions have been measured to 400/ J. A spin wave 
peak is located at w = 0.25 J and the shoulder like feature at oj ~ 0.5 J in 5/(q, w) 
(see Fig.||) is due to multi-spin-wave processes, the description of which is beyond 
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Fig. 4. Energy e(t) per spin for different order decomposition schemes for D = J: (solid line) 
predictor-corrector method; (dot-dashed line) 2nd-ordcr scheme; (dashed line) 4th-ordcr scheme; 
(dotted line) 8th-order method. The number of iterations performed are marked next to each line. 




Fig. 5. Dynamic structure factor of an isotropic Heisenberg ferromagnet for T = 0.8T C and 
q| = it/5 in the (100) direction on a simple cubic lattice (L = 10) obtained with the 2nd-order 
decomposition method for time step 8t = 0.04/ J. The longitudinal component is S;(q, u>) and the 
transverse component is St(q, u>). 
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the scope of this article. 
4. Summary 

We have described a set of algorithms which is based on Suzuki- Trotter decom- 
positions of exponential operators and compared their relative performance with 
each other as well as with a predictor-corrector method. The advantages of the 
predictor-corrector method are its versatility and its capability to conserve the mag- 
netization exactly. The decomposition of the lattice into sublattices, which is the 
basis for the decomposition method, depends on the range of the interactions, so 
that this approach is less general than the predictor-corrector method. Crystal 
field anisotropies leave the performance of the predictor-corrector method almost 
unaffected, whereas the decomposition method suffers from a drastic reduction in 
speed. 

The advantage of the decomposition method is its ability to handle large time 
steps and to conserve spin length exactly. In the absence of anisotropies it also con- 
serves the energy exactly and it maintains reversibility. For anisotropic Hamiltoni- 
ans energy conservation and reversibility can be obtained to a high accuracy using 
iterative schemes. Exact magnetization conservation, however, is lost. The time 
steps which can be used far exceed those used by the predictor-corrector method. 
In simple cases the 4th-order decomposition yields very accurate results even for 
time steps, which are an order of magnitude larger than typical time steps used for 
the predictor-corrector method. The 8th-order algorithm improves the conservation 
significantly but at the cost of greatly increased execution time. 
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